Two-Species Model System

This example is based on the system described in Exploratory Data Analysis Using Fisher Information (Mayer et al. 2007; Frieden and Gatenby 2007, 229–31). The two species Lotka-Volterra differential equations for the system are

\[ \frac{dx_1}{dt} = g_1 x_1 \left( 1- \frac{x_1}{k}\right) - \frac{l_{12} x_1 x_2}{1 + \beta x_1} \] \[ \frac{dx_2}{dt} = \frac{g_{21} x_1 x_2}{1 + \beta x_1} - m_2 x_2 \] The model parameters are \(g_1 = m_2 = 1\), \(l_{12} = g_{12} = 0.01\), \(k = 625\), and \(\beta = 0.005\). The initial conditions for the system were not provided in the original reference (Mayer et al. 2007). We found that \(x_1 = 277.7815\) and \(x_2 = 174.551\) provide reasonable results. The differential equations for the system are defined in R as shown below.

# Model parameters
parameters <- c(
  g1 = 1,
  m2 = 1,
  l12 = 0.01,
  g21 = 0.01,
  k = 625,
  B = 0.005
)

# Initial conditions
state <- c(
  x1 = 277.7815, 
  x2 = 174.551
)

# System differential equations (eq. 7.17 and 7.18)
deq <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dx1dt <- g1 * x1 * (1 - (x1 / k)) - (l12 * x1 * x2) / (1 + B * x1)
    dx2dt <- (g21 * x1 * x2) / (1 + B * x1) - m2 * x2
    list(c(dx1dt, dx2dt))
  })
}

We use the package deSolve to solve the system of differential equations. The original reference did not provide the time the system was observed (Mayer et al. 2007). We found that a complete cycle of the system corresponds to approximately 11.145 time units. The system of differential equations is solved and the system trajectory is plotted below.

# Vector of times
TT <-  11.145
times <- seq(0, TT, by = TT/1e3)

# Solve system differential equations
out <- ode(
  y = state,
  times
  = times,
  func = deq,
  parms = parameters,
  rtol = 1e-10,
  method = "ode45"
)

# Convert to data frame
sysSol <- data_frame(t = out[,1], x1 = out[,2], x2 = out[,3]) 
## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.
Phase space plot of system

Phase space plot of system

Change in predator and prey abundance over time

Change in predator and prey abundance over time

System State PDF & Dimensionality Reduction

The state of the system at any point in time, \(t\), is defined in two dimensions by the number of prey, \(x_1\), and the number of predators \(x_2\). Let \(\mathbf{x} = \left[x_1, x_2\right]\) represent the state of the system. The probability of observing the system in a particular state can be approximated based on a 2-dimensional histogram as shown below. Note that the system does not change at a constant rate (i.e., does not travel the path at a constant speed) and is therefore more likely to be observed in some states than others.

Normalized histogram showing probability of observing particular system states

Normalized histogram showing probability of observing particular system states

From top to bottom, distance traveled in phase space, speed tangential to system trajectory, acceleration tangential to system trajectory.

From top to bottom, distance traveled in phase space, speed tangential to system trajectory, acceleration tangential to system trajectory.

Comparison of normalized histogram of distance traveled (blue) to the derived probability density function (black)

Comparison of normalized histogram of distance traveled (blue) to the derived probability density function (black)

Fisher Information

The pdf of distance, \(p(s)\), was derived above. There are several equivalent equations for calculating shift-invariant FI. Some may offer numerical advantages over others. The general form is (Mayer et al. 2007, eq. 7.3b)

\[ I = \int \frac{ds}{p(s)}\left[\frac{dp(s)}{ds}\right]^2 \] The amplitude form is (Mayer et al. 2007, eq. 7.3c) \[ I = 4 \int ds\left[\frac{dq(s)}{ds}\right]^2 \]

A form specific to the pdf of distance traveled is derived as (Mayer et al. 2007, eq. 7.12) \[ I = \frac{1}{T} \int_0^T dt\left[\frac{s''^2}{s'^4}\right]^2 \] where integration is performed in the time domain. While all these formulations are equivalent, the last is most readily calculated when the differential equations for the system are known. Below all equations are shown to be equivalent. The FI for the distribution of distance is calculated to be approximately \(5.4 \times 10^{-5}\) which is consistent with published results (Mayer et al. 2007, 232).

Varying carrying capacity for the prey

Here we replicate the results of varying the prey carrying capacity, \(k\) (Mayer et al. 2007).

# Equation 7.3b
p <- sysSol_highO$p
s <- sysSol_highO$s
dp <- lead(p)-p
ds <- lead(s)-s
dpds <- dp/ds
ind <- 1:(length(s)-1)
FI_7.3b <- trapz(s[ind], (1/p[ind])*dpds[ind]^2)

# Equation 7.3c
q <- sqrt(sysSol_highO$p)
s <- sysSol_highO$s
dq <- lead(q)-q
ds <- lead(s)-s
dqds <- dq/ds
ind <- 1:(length(s)-1)
FI_7.3c <- 4*trapz(s[ind], dqds[ind]^2)

# Equation 7.12
t <- sysSol_highO$t
dsdt <- sysSol_highO$dsdt
d2sdt2 <- sysSol_highO$d2sdt2
ind <- 1:(length(s)-1)
FI_7.12 <- (1/TT)*trapz(t[ind], d2sdt2^2 / dsdt^4)

# Results
FI_7.3b
## [1] 5.370485e-05
FI_7.3c
## [1] 5.371751e-05
FI_7.12
## [1] 5.326461e-05
Phase space plot of system trajectories for different values of k

Phase space plot of system trajectories for different values of k

Plot of speed for different values of k

Plot of speed for different values of k

Plot distance pdf's for different values of k

Plot distance pdf’s for different values of k

Plot of Fisher information for different values of k

Plot of Fisher information for different values of k

Regime shift - Habitat destruction

Simulate a regime shift by abruptly decreasing the carrying capacity for the prey. The change in carrying capacity for the prey over time was modeled as

\[ k(t) = k_{1} - 0.5(k_1-k_2)(\tanh(\alpha(t - t^*)) + 1) \]

where \(k_1\) is the initial carrying capacity, \(k_2\) is the final carrying capacity, \(t^*\) is the time of the regime shift, and \(\alpha\) is a parameter that controls how abruptly the carrying capacity changes. The rate of change of the carrying capacity is given by

\[k'(t) = 0.5\alpha(k_1 - k_2)(\tanh(\alpha(t - t^*))^2 - 1) \]

Change in carrying capacity of prey over time for different values of alpha

Change in carrying capacity of prey over time for different values of alpha

Rate of change of carrying capacity of prey over time for different values of alpha

Rate of change of carrying capacity of prey over time for different values of alpha 

# Time of regime shift
tregime = 200

# Carrying capacity for regime 1 and regime 2
k1 <- 800
k2 <- 625

# Regime change rate
alpha <- 0.05
  
# Vector of times
TT <-  600
times <- seq(0, TT, by = TT / 2e3)

# Model parameters
parameters <- c(
  g1 = 1,
  m2 = 1,
  l12 = 0.01,
  g21 = 0.01,
  B = 0.005
)

# Initial conditions
state <- c(
  x1 = 61.9, 
  x2 = 35.3,
  s = 0,
  k = k1
)

# System differential equations including carrying capacity
deq_rs <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
  
    # Predator and prey rates of change
    dx1dt <- g1 * x1 * (1 - (x1 / k)) - (l12 * x1 * x2) / (1 + B * x1)
    dx2dt <- (g21 * x1 * x2) / (1 + B * x1) - m2 * x2
    
    # Rate of change of distance
    dsdt <- sqrt(dx1dt ^ 2 + dx2dt ^ 2)
    
    # Rate of change of carrying capacity
    dkdt <- 0.5*alpha*(k1 - k2)*(tanh(alpha*(t - tregime))^2 - 1)
  
    list(c(dx1dt, dx2dt, dsdt, dkdt))
  })
}

  
# Solve system differential equations (now with carrying capacity)
out_rs <- ode(
  y = state,
  times = times,
  func = deq_rs,
  parms = parameters,
  rtol = 1e-10,
  method = "ode45"
)
  
# Convert to data frame
sysSol <-
  data_frame(
  t = out_rs[, 1],
  x1 = out_rs[, 2],
  x2 = out_rs[, 3],
  s = out_rs[, 4],
  k = out_rs[, 5]
) %>%
  mutate(dsdt = (lead(s) - s) / (lead(t) - t),
  p = (1 / maxT) * (1 / dsdt)) %>%
  filter(!is.na(dsdt))
# Distance over which to move the window (in units time)
winspaceFactor <- 1

# Size of the window (in units time)
# winsizeVector <- seq(13.1, 11.0, by = -0.25)
winsizeVector <- c(13.061, 11.135)

winResults <- NULL
for(i in 1:length(winsizeVector)){
  
  # Define winsize
  winsize <- winsizeVector[i]
  
  for(j in 1:length(winspaceFactor)) {
    winspace <- winspaceFactor[j]*winsize
    
    # Start and stop points for windows
    winStart <- seq(min(times), max(times) - winsize, by = winspace)
    winStop <- winStart + winsize
    
    # Number of windows
    nWin <- length(winStart)
    
    FI_7.3b <- numeric(length(nWin))
    avgSpeed <- numeric(length(nWin))
    
    for (win in 1:nWin) {
      df <-
        sysSol %>%
        filter(t > winStart[win],
        t <= winStop[win]) %>%
        mutate(TT = max(t) - min(t),
        p = (1 / TT) * (1 / dsdt))
        
        # Fisher information
        p <- df$p
        s <- df$s
        dp <- lead(p) - p
        ds <- lead(s) - s
        dpds <- dp / ds
        ind <- 1:(length(s) - 1)
        FI_7.3b[win] <- trapz(s[ind], (1 / p[ind]) * dpds[ind] ^ 2)
        
        # Average speed
        avgSpeed[win] <- mean(df$dsdt)
    }
    
    tempResults <- data_frame(winStop = winStop,
    FI = FI_7.3b,
    avgSpeed = avgSpeed) %>%
    mutate(winsize = winsize,
           winspace = winspace,
           winspaceFactor = winspaceFactor[j])
    
    winResults <- rbind(winResults, tempResults)
    
  }
                            
}

winResults <- 
  winResults %>% 
  mutate(winsize = factor(winsize),
         winspace = factor(winspace))
Phase space plot of system trajectory with regime shift. Approximate location of regime shift is indicated by red dot.

Phase space plot of system trajectory with regime shift. Approximate location of regime shift is indicated by red dot.

Predator prey abundance with regime shift

Predator prey abundance with regime shift

Carrying capacity with regime shift

Carrying capacity with regime shift

Distance traveled in phase space with regime shift

Distance traveled in phase space with regime shift

Speed with regime shift

Speed with regime shift

Average speed over moving window with regime shift

Average speed over moving window with regime shift

Fisher Information over moving window with regime shift

Fisher Information over moving window with regime shift

Piecewise linear fit to distance traveled

Piecewise linear fit to distance traveled

Speed based on piecewise linear fit to distance traveled

Speed based on piecewise linear fit to distance traveled

References

Frieden, B. Roy, and Robert Gatenby, eds. 2007. Exploratory Data Analysis Using Fisher Information. Springer. http://www.springer.com/us/book/9781846285066.

Mayer, Audrey L., Christopher W. Pawlowski, Brian D. Fath, and Heriberto Cabezas. 2007. “Applications of Fisher Information to the Management of Sustainable Environmental Systems.” In Exploratory Data Analysis Using Fisher Information.